clear frames
use "data/ready4analysis.dta",clear

xtset pidp year
							  
****************************************************************************************
**** FIGURE 3
****************************************************************************************

gen age_mc=age-40
recode birthy (1935/1943=1) (1944/1952=2) (1953/1960=3) (1961/1968=4) (1969/1976=5) (nonm=.), gen(gen5)

* Covariates for main model and model in supporting information
local covars_main "i.maxeduc" // Main model

foreach dv in z { // Change to "a b c d e f" for individual items
	foreach type in lifetime  { // change to panel for panel stable class
	foreach mod in main { // Optional for alternative covariates
	foreach t in 75  { // Change to 90 for more strict class classification
	local model "gen3" // Change to gen5 to use 5 generation model 
		
	if "`model'"=="gen3" {	
	* USING MAIN COHORT SPECIFICATION 	
		eststo clear
		eststo fig2: mixed red`dv' i.egp_`type'`t' i.egp_`type'`t'#c.age_mc c.age_mc##c.age_mc i.polgen i.polgen#c.age_mc i.female `covars_`mod'' ///
								  || _all:R.year || pidp:age_mc if inrange(birthy,1935,1976) & ///
								  inrange(N_red,5,7) , cov(unstruct)  stddev
			
			matrix N_g = e(N_g)
			local groups = N_g[1,1]
			estadd scalar groups = `groups'

		
		esttab fig2 using "tables/figure2_`dv'_egp_`type'_`mod'_`t'.rtf", replace b(3) se(3) transform(ln*: exp(@) exp(@) at*: tanh(@) (1-tanh(@)^2))         ///
			 drop(1.egp_`type'* 1.polge* 1.maxe* 0.female) eqlabels("" "sd(age)" "sd(_cons)" "cov(age,_cons)" "sd(Residual)", none) ///
			 varlabels(,elist(weight:_cons "{break}{hline @width}"))   ///
			 varwidth(30) stats(N groups aic bic, labels("N" "Groups" "AIC" "BIC"))

		margins,at(polgen=3 age_mc=(-24(2)6) egp_`type'`t'=(1 2 3 4 5 120)) ///
				at(polgen=2 age_mc=(-7(2)20) egp_`type'`t'=(1 2 3 4 5 120)) ///
				at(polgen=1 age_mc=(7(2)25)  egp_`type'`t'=(1 2 3 4 5 120)) post 

		cap graph drop `dv'egp
		marginsplot, x(age_mc) by(egp_`type'`t') byopt(r(1)) name("`dv'egp")

		graph export "temp/fig2_`dv'_egp_`mod'_`type'`t'.png", as(png) name("`dv'egp") replace

		// Below saves results for R
		qui parmest, saving("temp/fig2.dta",replace) 

		frame create new
		frame change new
			
			use "temp/fig2.dta", clear
				
			matrix varnames = e(at)
			svmat varnames
			
			gen gen=varnames8+2*varnames9+3*varnames10
			gen class=varnames1+2*varnames2+3*varnames3+4*varnames4+5*varnames5+6*varnames6
			rename varnames7 age
				replace age=age+40

			save "figs/fig2_`dv'_egp_`mod'_`type'`t'.dta", replace

		frame change default
		frame drop new
				
		save "temp/fig2_data_`dv'_egp_`mod'_`type'`t'.dta", replace
	}
	else {
		
	* USING ALTERNATIVE COHORT SPECIFICATION 
	eststo clear		 
	eststo fig2: mixed redz i.egp_`type'`t' i.egp_`type'`t'#c.age_mc c.age_mc##c.age_mc i.gen5 i.gen5#c.age_mc i.female `covars_`mod'' ///
								  || _all:R.year || pidp:age_mc if inrange(birthy,1935,1976) & ///
								  inrange(N_red,5,7) , cov(unstruct) stddev 
	  
		matrix N_g = e(N_g)
		*local groups = N_g[1,1]
		*estadd scalar groups = `groups'
		estadd scalar groups = N_g[1,1]
		
			
		esttab fig2 using "tables/figure2_`dv'_egp_`type'_`mod'_`t'_`model'.rtf", replace b(3) se(3) transform(ln*: exp(@) exp(@) at*: tanh(@) (1-tanh(@)^2))         ///
		 drop(1.egp_`type'* 1.gen* 1.maxe* 0.female) eqlabels("" "sd(age)" "sd(_cons)" "cov(age,_cons)" "sd(Residual)", none) ///
		 varlabels(,elist(weight:_cons "{break}{hline @width}"))   ///
		 varwidth(30) stats(N groups aic bic, labels("N" "Groups" "AIC" "BIC"))

		margins,at(gen5=5 age_mc=(-24(2)-6) egp_`type'`t'=(1 2 3 4 5 120)) ///
			at(gen5=4 age_mc=(-14(2)2) egp_`type'`t'=(1 2 3 4 5 120)) ///
			at(gen5=3 age_mc=(-6(2)10) egp_`type'`t'=(1 2 3 4 5 120)) ///
			at(gen5=2 age_mc=(2(2)18) egp_`type'`t'=(1 2 3 4 5 120)) ///
			at(gen5=1 age_mc=(10(2)26)  egp_`type'`t'=(1 2 3 4 5 120)) post 

		cap graph drop egp
		marginsplot, x(age_mc) by(egp_`type'`t') byopt(r(1)) name("egp")

		graph export "temp/fig2_egp_`type'`t'_`mod'_`model'.png", as(png) name("egp") replace

		qui parmest, saving("temp/fig2.dta",replace) 

		frame create new
		frame change new
			
			use "temp/fig2.dta", clear
				
			matrix varnames = e(at)
			svmat varnames
			
			gen gen=varnames8+2*varnames9+3*varnames10+4*varnames11+5*varnames12
			gen class=varnames1+2*varnames2+3*varnames3+4*varnames4+5*varnames5+6*varnames6
			rename varnames7 age
			replace age=age+40

			save "figs/fig2_z_egp_`type'`t'_`mod'_`model'.dta", replace

		frame change default
		frame drop new

		}
		}
	}	
	}
	}

		